function [TCP, AA] = calculate_TCP(param, delta_t, free, t_f1, t_f2, param_id, sigma, n, LQL, vd, D, t_rad, t_treat_c4, use_Markov, AA)

% INFO: Function to calculate TCP

%% OUT:

    % TCP: (real) TCP
    % AA: (matrix) parameters used for each trial
    
%% IN:

    % param: (vec) parameters
    % delta_t: (real) time step
    % free: (vec) controls free growth and start of treatment condition
    %    free(1): (int) 1 for free behavior until an specific time or volume, 0 otherwise 
    %    free(2): (int) 1 for time, 0 for volume
    %    free(3): (real) initial volume or time
    % t_f1: (real) maximum time allowed to reach the initial state
    % t_f2: (real) simulation time after initial state is reached
    % param_id: (vec) variation parameters index
    % sigma: (real) parameters variation
    % n: (int) number of experiments to calculate the TCP
    % LQL: =1 if LQL, otherwise (modified) LQ (with PARAM(33)) 
    % vd: =1 if vascular death (D > 15)
    % D: (vec) RT fractions
    % t_rad: (vec) RT treatment days
    % t_treat_c4: (vec) IT treatment days
    % use_Markov: =1 if Markov, otherwise deterministic
    % AA: (matrix) parameters to each trial or 0 to generate new ones

%% Algorithm

control = 0;
val_cases = 0;
flag = 1; %The parameters for each trial are given
if AA == 0 
    flag = 0; %Generate new parameters
    AA = zeros(n, length(param_id));
end

%Gaussian perturbations
mu = param(param_id);
sigma = sigma .* mu;
                    
ii = 1;
trials = 0;

while val_cases < n
    trials = trials + 1;
    if flag == 1
        param(param_id) = AA(ii,:);
    else
        A = normrnd(mu, sigma);
        A(A < 0) = 0; 
        param(param_id) = A;
    end
    
    if t_treat_c4 == 0
        param(25) = 0;
    end

    [ ~, t_eq, ~, ~, sol] = radioimmuno_response_model(param, delta_t, free, t_f1, t_f2, D, t_rad, t_treat_c4, 0, LQL, vd, use_Markov);

    if t_eq == -1 % Non-valid simulation
    else
        if flag == 0
            AA(ii, :) = A;
        end
        ii = ii + 1;
        sol_eq = sol(round(t_eq/delta_t):round(t_eq/delta_t + t_f2/delta_t));

        val_cases = val_cases + 1;
        
        if (sol_eq(length(sol_eq)) == 0) % Using Markov
            control = control + 1;
        end
    end
    
    if trials > 2*n
        TCP = 0;
        return
    end
end %while

%Calculate TCP
TCP = max(0, control/val_cases);

end %function